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We present the first numerical approximation of the scalar Schwarzschild Green function in the 
time domain, which reveals several universal features of wave propagation in black hole spacetimes. 
We demonstrate the trapping of energy near the photon sphere and confirm its exponential decay. 
The trapped wavefront passes through caustics resulting in echoes that propagate to infinity. The 
arrival times and the decay rate of these caustic echoes are consistent with propagation along 
null geodesies and the large I limit of quasinormal modes. We show that the fourfold singularity 
structure of the retarded Green function is due to the well-known action of a Hilbert transform 
on the trapped wavefront at caustics. A twofold cycle is obtained for degenerate source-observer 
configurations along the caustic line, where the energy amplification increases with an inverse power 
of the scale of the source. Finally, we discuss the tail piece of the solution due to propagation within 
the light cone, up to and including null infinity, and argue that, even with ideal instruments, only 
a finite number of echoes can be observed. Putting these pieces together, we provide a heuristic 
expression that approximates the Green function with a few free parameters. Accurate calculations 
and approximations of the Green function are the most general way of solving for wave propagation 
in curved spacetimes and should be useful in a variety of studies such as the computation of the 
self-force on a particle. 



I. INTRODUCTION 

A general technique to study the full response of a 
black hole to generic perturbations is to explore the 
Green function. While the Green function (aptly called 
the propagator) is simple in flat spacetime, it has a rich 
structure in curved spacetimes. 

Recently, Ori discovered that the Green function in 
Schwarzschild spacetime displays a fourfold periodic 
structure [1]. This fourfold structure can be understood 
in terms of trapped null geodesies, similar to the clas- 
sical interpretation of quasinormal mode ringing as en- 
ergy leakage from the photon sphere [5J [3] . At each half 
revolution along the trapped orbit, the trapped wave- 
front passes a caustic undergoing a transformation with 
a fourfold cycle. Ori investigated this cycle using a quasi- 
normal mode expansion method. He also performed an 
acoustic experiment providing evidence for the structure 
of trapped signals through caustics. 

The fourfold cycle has been formally analyzed in Gen- 
eral Relativity for the first time by Casals, Dolan, Ot- 
tewill, and Warded in |H [S] in the example of Nariai 
spacetimes. Recently, Casals and Nolan developed a 
Kirchhoff 's integral representation of the Green function 
on Plebahski-Hacyan spacetimes, confirming the appear- 
ance of the fourfold singularity structure in these black- 
hole models [6]. Going beyond toy models, Dolan and 
Ottewill discussed the Green function in Schwarzschild 
spacetime using large t quasinormal mode sums [7] based 
on their expansion method [5j. The essential singular- 
ity structure of the retarded Green function in arbitrary 
spacetimes has been discussed by Harte and Drivas [S]. 
They approximate the region of any spacetime near a null 
geodesic by a pp-wave spacetime in the Penrose limit. 
Their analysis implies, in particular, that the fourfold 
structure should be valid also in Kerr spacetimes. 



Motivated by these theoretical studies and by Ori's 
acoustic experiment, we perform a numerical experiment 
and analyze in detail the response of a black hole to com- 
pactly supported scalar perturbations. We describe the 
problem and our numerical setup in Sec. [XXJ. The main 
part of this paper (Sec. Ill I contains our results orga- 



nized by scales, from the geometrical optics limit (short 
wavelength) to the curvature scale (long wavelength). 
After a qualitative description of the evolution of com- 



pactly supported wave packets in Sec. Ill A we discuss 
the propagation of the signal along null geodesies and the 
subsequent caustic echoes as measured by an observer 



at infinity (Sec. Ill B I . Wave propagation through caus- 



tics and the fourfold structure of the Green function is 
studied in Sec. IIII CI We show that the fourfold struc- 
ture is due to the action of a Hilbert transform at each 
passage through a caustic. In Sec. |III D| we present a 
twofold cycle for degenerate source-observer configura- 
tions and find an inverse power relation between energy 
magnification at caustics and the source scale. We ar- 
gue that the trapping of energy near the black hole is 
the universal cyclic feature manifested in the Green func- 
tion, not the number of the cycle (e.g., four or two). In 
Sec. |III E| we analyze the tail of the signal due to scat- 
tering off the background curvature, which dominates at 
late times implying that only a finite number of echoes 
can be observed even with ideal instruments. Putting 
these pieces together in Sec. |IIIF| we provide a heuristic 
expression for the Green function with a few free param- 
eters. We present our findings, discuss the limitations 
of our method, and pose problems for future research in 
Sec. |IV| The Appendixes complement the main text. The 
topics include hyperboloidal compactification for numer- 
ical simulations (Appendix A]), the Schwarzschild Green 
function in the geometrical optics limit (Appendix [B| , 
and the Hilbert transform at caustics (Appendix [C| . 
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II. THE SETUP 

A. The problem description 

The simplest model equation to describe essential fea- 
tures of wave propagation is the scalar wave equation 



D<f>(x) = S(x) , 



(1) 



for a scalar field (f> with source S that depends on space- 
time coordinates x. The wave operator, □, is written 
with respect to a background Schwarzschild metric. We 
set the source to 



S(x) 



1 



(27TCT 2 ) 



2 eX P 



<^„(y - x'»)(x v - X ,v ) 
2a 2 



, (2) 



where 5^ is the four-dimensional Kronecker delta and 
a sets the scale for the width of the perturbation. The 
Gaussian source to the scalar wave equation is chosen to 
approximate the delta distribution because we are inter- 
ested in the evolution of wavefronts and in the numerical 
approximation of the Green function. The Schwarzschild 
Green function to the scalar wave equation satisfies 



-in 



DG(x,x') = — 5 i (x-x'), 



(3) 



where g is the metric determinant, and 8 is the four- 
dimensional Dirac delta distribution. The solution to 
the wave equation with a narrow Gaussian source ^ 
provides an approximation to the Schwarzschild Green 
function with a delta distribution source up to an overall 
factor. 



B. The numerical setup 

We numerically solve the scalar wave equation |lj with 
source ^ using the Spectral Einstein Code SpEC [ID] . 
SpEC is a spectral element code for solving elliptic and 
hyperbolic partial differential equations with particular 
focus on the Einstein equations for the simulation of com- 
pact binaries. 

The time evolution of hyperbolic equations in SpEC 
uses the method of lines. A spectral expansion of the 
unknown in space is performed in elements that commu- 
nicate with each other along internal boundaries through 
the exchange of characteristics via penalty terms. The 
discretized unknown is then evolved as a coupled system 
of ordinary differential equations in time with adaptive 
time-stepping using the Dormand-Prince method. 

The topology of the grid in SpEC can be chosen de- 
pending on the particular problem. We solve the scalar 
wave equation on a Schwarzschild background with in- 
ner boundary at the event horizon and outer boundary 
at null infinity. Therefore, our grid consists of concentric 
spherical shells that span the domain between the hori- 
zon and infinity. We use Chebyshev polynomials with 



Gauss-Lobatto collocation points in the radial direction 
and a spherical harmonic expansion in the angular direc- 
tion. The topology of the grid and the initial pulse ^ 
are depicted in Fig. [I] 




FIG. 1: Topology of the numerical grid and the initial pulse 
perturbation on the equatorial plane. We only show 5 concen- 
tric shells of width 2 spanning the radial coordinate domain 
p £ [1.99, 12] with an expansion into spherical harmonics of 
order 47. The simulations presented in Sec. |III| use 40 such 
shells with a spherical harmonic expansion of order 87. 



We employ hyperboloidal scri-fixing [11] to compute 
the unbounded domain solution and the signal at infinity 
as measured by idealized observers. The coordinates are 
constructed using the hyperboloidal layer method [H] 
based on the ingoing Eddington-Finkelstein coordinates 
in the interior (Appendix [A|. Setting the Schwarzschild 
mass M, which provides the scale for the coordinates, 
to unity, the coordinate domain in the radial direction 
is p e [1.99, 12] where p = 2 corresponds to the event 
horizon and p — 12 corresponds to null infinity. The 
interface to the hyperboloidal layer is at p = 8. The 
hyperboloidal layer consists of the two outermost shells 
on the grid in Fig. [T] 

The simulations presented in the results section 
(Sec. Ill) employ 40 concentric shells with 7 collocation 
points each and a spherical harmonic expansion of de- 
gree 87. For Fig. [I] we use a sparse grid for clarity 
of visualization. The Gaussian pulse ^ is based at 
x in _ {t' )X ',y',z'} = {2,6,0,0}. Our results are robust 
under variation of parameter values. 



We set the standard deviation as a = 0.2 unless oth- 
erwise stated, thus giving a 10:1 ratio of the radius of 
the event horizon to the scale of the source. We would 
obtain the Green function in the limit of vanishing width 
of the Gaussian and infinite numerical resolution. Our 
numerical approximation can be improved using smaller 
a and higher resolution. 
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FIG. 2: Top: The evolution of the scalar field generated by 
a compact source at 6M on the positive x axis as measured 
by an observer at future null infinity on the positive z axis. 
Bottom: The absolute value of the same field on a log scale, 
showing the detailed structure in the amplitude. 



III. RESULTS 
A. Qualitative description 

It is well-known that a generic scalar field in 
Schwarzschild spacetime dissipates to infinity in an expo- 
nentially decaying ringing and a polynomially decaying 
tail. The fast dissipation is confirmed in the top panel 
of Fig. [2j which shows the time evolution of the field as 
measured by an observer on the positive z axis at null in- 
finity. Following the initial pulse, we see more structure 
that becomes visible on a log scale in the bottom panel. 
There seems to be a recurring signal caused by the ini- 
tial perturbation that disappears at late times when the 
polynomial decay dominates. 

The dynamics behind Fig. [2] can be seen more clearly 
in three-dimensional plots. In Figs. [3][6]we show snap- 
shots from the evolution of <p on the equatorial plane. 
We set the width of the Gaussian source to a = 0.3M for 
visualization. 

The curvature of spacetime near the black hole causes 
the wavefront triggered by the initial pulse perturbation 
in Fig. [T] to bend and to focus forming a half-cardioid 
shape typical of caustic formation (Fig. [3]) . At the cusp of 
the cardioid, by the antipodal point to the initial pertur- 
bation near the black hole, the wavefront passes through 
itself thus forming a caustic. The caustic moves out to 
infinity along the negative x axis and leaves an echo in its 
wake. This secondary wavefront again forms a cardioid 
shape with the cusp on the other side of the black hole 
(Fig. [2]). This process repeats itself resulting in echoes 
that reach observers in regular time intervals (Figs.pland 

§■ 

The figures suggest that the echoes are due to par- 
tial trapping of the initial energy at the photon sphere 
at areal coordinate 3Af. Part of the wavefront from the 
initial perturbation dissipates directly to infinity, part of 




FIG. 3: Formation of the first caustic echo from the initial 
pulse in Fig. [I] We plot a generic bird's eye view onto the 
equatorial plane. A movie of the dynamical evolution, includ- 
ing a top to bottom view, can be seen online [13]. The color 
and height of the surface is determined by the amplitude of 
the scalar field <j>. The wavefront of the direct signal has a 
Gaussian profile and forms a cardioid shape around the event 
horizon, which is typical of caustic formation. The amplitude 
amplification is clearly visible. The first caustic echo forms 
in the wake of the caustic and is seen as a wavefront in the 
vicinity of the horizon. It generates a second caustic echo 
plotted in Fig. [4] 




FIG. 4: Formation of the second caustic echo. The profile of 
the signal from the wavefront of the first echo, seen at the 
boundary, is not a Gaussian but its Hilbert transformation 



see Sec. Ill C I. The wavefront of the second caustic echo near 



the horizon is strictly negative and has a negative Gaussian 
profile (see boundary profile in Fig. [5|. 



it falls into the black hole, and part is trapped near the 
horizon forming so-called "surface waves" propagating on 
the photon sphere [3J |3J Q31 [TS] . The trapped signal goes 
through two caustics upon each full revolution and leaks 
out energy to infinity in the form of propagating wave- 
fronts that we call caustic echoes. 

The profiles of the caustic echoes follow a certain pat- 
tern. The direct signal has the profile of a Gaussian in 
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FIG. 5: Formation of the third caustic echo. The second 
caustic echo has a negative Gaussian profile as seen at the 
boundary. The profile of the third echo is clearer in Fig. [6] 




FIG. 6: Formation of the fourth caustic echo. The profile of 
the third caustic echo is the Hilbert transform of the nega- 
tive Gaussian, which is the negative of the profile of the first 
caustic echo seen in Fig. [4] The fourth caustic echo emanat- 
ing from the horizon has a positive Gaussian profile like the 
direct signal, thus completing the fourfold cyle. 



accordance with the source (Figs. [T] and |3j. The shape of 
the first caustic echo (see the boundary profile in Fig. [4]) 
is the Hilbert transform of the direct Gaussian signal, 
which we refer to as a Dawsonian, and is discussed in 
more detail in Sec. |III C| The second echo is a negative 
Gaussian (Fig. [5]) and finally the fourth echo is a negative 
Dawsonian (Fig. |6j. The wavefront emanating from the 
wake of the fourth echo has a Gaussian profile just like 
the initial signal, thereby completing the fourfold cycle 
(Fig. §• 

This is the demonstration of the fourfold structure first 
discovered in General Relativity by Ori [T] and analyzed 
by Casals et al. [4J [5] . A visualization of the time evolu- 
tion can be found on the internet |13j . In the following, 
we present a detailed quantitative understanding of the 
data presented in Fig. [2] We present analysis for the 



measurements of an observer at future null infinity on 
the z axis, unless stated otherwise. 



B. Arrival and decay of caustic echoes 

The first feature seen by the observer in Fig. [2] is a nar- 
row Gaussian pulse that we refer to as the direct signal, 
which is the part of the perturbation that propagates di- 
rectly to the observer. We set the observer's clock to zero 
by the arrival time of the direct signal's maximum. The 
arrival times of the subsequent caustic echoes are plotted 
in Fig. 
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FIG. 7: Time differences between caustic echoes indicating 
agreement with theoretical prediction based on null geodesic 
revolution times around the photon sphere at 3M. Our mea- 
surement of the arrival times for differently shaped echoes 
leads to an artificial, systematic deviation from the theoreti- 
cal prediction (ttv27M, dashed line), which averages out (left 
panel). When similarly shaped echoes are compared, however, 
the measured time differences (between even echoes in the 
middle panel and between odd echoes in the right panel) agree 
accurately with the revolution time of null geodesies around 
the photon sphere given by Tf u u = 2ny/27M w 32.65A/. 

Even and odd echoes, that is, echoes resulting from 
the wavefront passing through an even or odd num- 
ber of caustics, have different profiles (see Figs. 2][6 and 



Sec. Ill C ). Therefore, it is appropriate to consider the ar- 



rival times of even and odd echoes separately. We extract 
the times between even echoes by fitting a Gaussian and 
reading off the time at which the peak amplitude arrives 
at the observer (middle panel of Fig. [7]) . We extract the 
times between odd echoes by fitting their time derivative 
to a Gaussian and finding the time at which the peak of 
the derivative is maximal (right panel of Fig. [7]) . 

The arrival times of caustic echoes approach a con- 
stant value. This observation is explained through the 
generation mechanism behind the echoes, namely, the 
revolution of trapped energy at the photon sphere along 
null geodesies. The quantitative prediction of the ar- 
rival times requires only the study of high-frequency 
wave propagation for which it suffices to employ the ge- 
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ometrical optics approximation. An analysis of solutions 
to the scalar wave equation in Schwarzschild spacetime 
with localized energy has been performed by Stewart 
[IB] . Stewart shows that the peak amplitude propagates 
along null geodesies, which orbit the photon sphere with 
a period of T M1 = 2nVffM « 32.6484M. The half- 
period for the formation of each successive caustic is then 
T haJf = ttV27M « 16.3242M. These values are plotted 
as dashed lines in Fig. [7] and indicate good agreement of 
the measured arrival times with the theoretical predic- 
tion. The best-fit value for the successive echoes shown 
in the left panel of Fig.[7]agrees with Thaif to 4 significant 
digits. 

The geometrical optics approximation is a short-wave 
(or high-frequency) approximation, and should therefore 
agree with the large I limit of quasinormal modes. In this 
limit, the quasinormal mode frequencies can be expressed 
analytically as (TT1 [18] 

Mu%r i ~J=(t-i( n + 1/2)). (4) 

The real part of the frequency (divided by £) corresponds 
to the full orbital frequency 27r/Tf u n of null rays on the 
unstable photon orbit. Thus, the time between successive 
even or odd caustic echoes is related to the real part of 
quasinormal mode frequencies in the large £ limit, which 
in turn are related to properties of null geodesies at or 
near the unstable photon orbit. 

From Q it is clear that the amplitude of caustic echoes 
should decay exponentially with a decay rate given by 
the imaginary part of the large I frequency for the lowest 
overtone (hence longest decay time), n — 0. The same 
decay rate can be calculated from the geometrical op- 
tics limit [TB] and is given by the Lyapunov exponent of 
unstable photon orbits [31 [05] . 




t/M 

FIG. 8: Amplitude decay of the field in Fig. [2] The dashed 
line shows the exponential decay implied by the imaginary 
part of the n = overtone of the quasinormal mode in the 
large £ limit Q, or equivalently the Lyapunov exponent A of 
the null geodesies near the unstable photon orbit. 

Figure [2] indicates that the amplitude of each caustic 
echo decays exponentially with time until the backseat - 



ter off the background curvature dominates. Zooming 
into the interval before the tail-dominated domain, Fig. [8] 
shows visual agreement between the data and the the- 
oretically predicted decay rate of A = l/(2\/27)M _1 « 
0.096225A/ -1 . For comparison, the decay rate associated 
with the peak amplitudes of the even (odd) caustic echoes 
obtained by fitting to Ae~ x ™™ 1 gj ves \ num = 0.103M -1 
(0.097M- 1 ). 

The arrival and decay times of caustic echoes ap- 
proach their predicted values from the geometrical optics 
approximation (or equivalently the large £ quasinormal 
mode limit) very fast. The first few echoes following the 
direct signal, however, deviate from these limiting values 
considerably (Fig. [7|) . This deviation is due to the propa- 
gation of the first caustic echo along a null geodesic that 
is farther away from the photon sphere than the subse- 
quent echoes. 




FIG. 9: Cartoon depicting several null geodesies connecting 
the source (yellow star) to the observer (stick figure) for a 
generic relative angle, 7 £ [0, n). The black hole is represented 
by a solid black disk and its unstable photon orbit by a dotted 
circle. Left: The null geodesic paths corresponding to the 
first two odd caustic echoes. Right: The null geodesic paths 
corresponding to the first three signals. 

A schematic diagram of the propagation of the first 
three signals is given in Fig. [9] We label the echoes 
by the order in time in which they appear, or equiva- 
lently by the number of caustics N that the wavefront 
has passed through before leaving the photon sphere. 
Consider the the first pulse after the direct signal that 
reaches the observer, the N = 1 echo. The pulse travels 
along the null geodesic connecting the source to the ob- 
server wrapping partially around the far side of the black 
hole. The angle swept through by this null geodesic is 
Aipi = 3ir/2. The null geodesic corresponding to the 
N = 3 echo makes a full orbit around the black hole so 
that the angle swept through is Aip^ — 7tt/2. For any rel- 
ative angle between the source and observer, 7 6 [0, n], it 
is clear from Fig. [9] that the time between successive even 
or odd caustic echoes approaches the period for a null 
geodesic to make a full orbit around the photon sphere; 
this is a universal property of the black hole independent 
of the source-observer configuration. 

For generic source-observer configurations, the time 
between caustic echoes depends on the angle 7 between 
the source and observer. The times between successive 
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caustic echoes are given by 



7 



T, 



full 



N even 



AT N ( 7 ) = 



7T — 7 



(5) 



Tsui , iVodd 



The cases 7 = 0, 7r are degenerate because AT/v = for 
either an even or odd echo. These observers see a twofold 
cycle as opposed to a fourfold cycle, and the peak ampli- 
tude of the signal is amplified as seen in Figs. [3}|6] The 
amplification cannot be treated in the geometrical optics 
approximation because the amplitude at caustic points 
is unbounded in this limit. We discuss the twofold cycle 
and the energy amplification for degenerate observers in 
Sec. HHPl 



C. Profiles of caustic echoes 

The profiles of caustic echoes transform by propagating 
through caustics, which leads to a fourfold cycle [UIHE]. 
In this section, we analyze in detail the structure of the 
pulse profiles. 



Figure 10 shows the profiles for the same source- 
observer configuration as in the previous section. The 
first pulse in the upper left panel is the direct signal and 
the subsequent pulses are the caustic echoes. The four- 
fold cycle presented in each row repeats itself until the 
signal is dominated by the tail (third row of Fig. 10 1. 



More caustic echoes are visible in the local decay rate 



(Fig. 151 



Using the geometrical optics approximation in Ap- 
pendix[Cj we show that passing through a caustic induces 
a phase shift of the wavefront in the frequency domain 
by — 7r/2 radians and is equivalent to a Hilbert transform 



( C4 ) . The phase change of light traveling through a caus- 
tic is a well-known effect in optics, called the Gouy phase 
shift, and was discovered in 1890 [5D] (for a more recent 
discussion see Chapter 9.5 of [H]). The fourfold struc- 
ture of the Green function is then a simple consequence 
of trapped null geodesies passing through caustics with 
each passage causing a phase shift by — tt/2, completing 
a full cycle of 2ir in four echoes. 

The Green function in the geometrical optics approxi- 
mation is given by ( C7 1 



G(x,x')~ £ |4)(x,x')|0(t-t / -T A r_i(x s x')) 

Re{H N [6(t-t' -7Mx,x'))]}. (6) 



JV=0 



Here, Tjv(x, x') is the coordinate time along a null 
geodesic connecting the spatial points x and x' and pass- 
ing through N caustics, Hjy denotes N applications of 
the Hilbert transform (thus phase shifting its argument 
by —Nn/2), the step function ensures causality (with 
T_i(x,x') = 0), and A (x, x') is the leading order field 
amplitude in the geometrical optics limit (see Appendix 
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FIG. 10: The first 12 pulses (direct signal plus 11 caustic 
echoes) seen by an observer on the positive z axis at future 
null infinity. The profiles of the caustic echoes are repeated 
in a fourfold cycle as is clearly seen in the first two rows. In 
the last row, the tail effect from backscattering dominates the 
signal. 



[C] for the derivation) . The field is obtained by the con- 
volution of the Green function with the source S(x) 



~ 00 
cj>{x)= G{x,x')S{x')^Y,^ x ^ 



(7) 



where f , = J d 4 x' —g(x') with g the determinant of 
the spacetime metric, and 4>n{x) is the field evaluated 
in the geometrical optics approximation associated with 
the passage of the wavefront through ./V caustics. 

Equation Q gives us a quantitative prescription for 
the profiles of the wavefronts after propagation through 
N caustics. For N — 0, the Hilbert transform Hq returns 
its argument. Hence the direct signal obtained from the 
convolution of the approximate Green function (|6| with 
the Gaussian source |2| is proportional to a Gaussian as 
confirmed in the upper left panel of Fig. [To} For N = 1 
the Hilbert transform Hi applies a phase shift by — tt/2 
radians to the signal entering the caustic. For the delta 
distribution in ^ this is given by (C3|. Interestingly, 
( |C3[ ) has a different singularity structure than the delta 
distribution and has support in the interior of the forward 
light cone of x'. Hence, this contribution adds to any tail 
piece from backscattering once the wavefront has passed 
through an odd caustic. The caustic thus smears out the 
(sharp) delta distribution that gets refocused to a delta 
distribution when passing through the next caustic. 
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For the first caustic echo we get from Q 



, le-3 
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xRe 



1 [ \ A ( 'M ( ( x ' ~ x o)' 
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du ( iix , . . 

2tt 6XP V "2" ~ l ~ X 

.t-T^x.x') / 

x / lit' exp I iwi' 

J —00 \ 



2a 2 



(8) 



Integrating over t' (assuming t — Ti(x, x') > £o + er so 
that the upper integration limit is replaced by 00) and 
oj, and assuming that the spatial part of the Gaussian is 
sufficiently narrow to be replaced by a delta distribution 
yields 



h(x) Pa -|A (x,x )|^(i-to-T 1 i(x,xo),a), (9) 



where 



D(y,a) = e- y2 /( 2CT2 )erfi 



(10) 



We will refer to ( 10 ) 



is related to Dawson's integral [22] 
as the "Dawsonian" below. 

The upper left panel of Fig. [TT] shows the first caustic 
echo (e.g., the second pulse in the upper left panel of 
Fig. [l0| together with two different fits to a Dawsonian 
using (I9J). The black line is the numerical data and the 
red line is a fit to aD(t—b, c)+d with a = — 2.8x 10 -3 , b — 
38.64M, c = 0.28M, and d = 2.4x 10~ 5 . Note that the 
best-fit value for c = 0.28M differs substantially from the 
original width of the Gaussian source of a — 0.2M. This 
free fit is remarkably good over a large interval around 

r 1 (x,x ) = 6. 

The black dashed line in Fig. [TT] shows the approxima- 
tion |9]) with <t = 0.2M (i.e., the width of the Gaussian 
source) and Ti(x,Xo) fixed to the value extracted from 
the numerical simulation; only the amplitude is fitted for. 
Despite the crudity of the approximations, the restricted 
fit captures qualitative features of the N = 1 echo rather 
well. The agreement indicates that the Green function 
in the geometrical optics approximation is suitable to de- 
scribe the propagation of wavefronts through caustics in 
terms of Hilbert transformations. The remaining panels 
in Fig. 11 for the next three caustic echoes show similar 



results. 

Our discussion so far covers wavefront propagation 
through caustics for relatively narrow pulses providing 
an approximation to the Green function. We confirmed 
the action of a Hilbert transform at caustics using other 
waveforms as well. 

The top panels in Fig. [12] show a fourfold cycle trig- 
gered by an initial square profile, that is, a pulse localized 
in space but persisting for a longer time scale than the 
time scale of the black hole. In the lower panel of Fig. |T2"] 
we plot the direct signal (black curve) and compare its 
Hilbert transform (blue curve) to the rescaled and shifted 




le-4 



Data 
Free fit 
Restricted fit 







57 



62 
le-5 



68 73 
t/M 



85 90 

t/M 



1 

-1-2 
■3 
■4 
■5 
2.0 
1.5 
1.0 
0.5 
0.0 
0.5 



95 



FIG. 11: The first four caustic echoes in the numerical solu- 
tion (black) along with two fits to a Dawsonian using the ge- 
ometrical optics approximation (red and dashed black). The 
numerical solution is fit to Eq. |9| with the amplitude, time 
of arrival, and width as fitting parameters. The fit is so well 
that the red and the black lines are hardly distinguishable. 
The black dashed line is a restricted fit where only the am- 
plitude is fitted for and the time of arrival and the width are 
set by the numerical solution and the Gaussian source width, 
respectively. The quality of the fits indicate that the geomet- 
rical optics approximation captures the profiles of the caustic 
echoes accurately via the Hilbert transform. 



profile of the N — 1 caustic echo from the numerical 
simulation (red curve). The agreement between the blue 
and the red curves implies that the wavefronts are indeed 
Hilbert transformed through caustics. 



D. The twofold cycle and the caustic line 

Equation ^ indicates that for angles 7 = 0, it be- 
tween the source and the observer, the structure of the 
measured signal is degenerate. This is the caustic line 
where caustic echoes merge (the x axis in our simula- 
tions). Even and odd echoes arrive at the same time at 
the caustic line, therefore we obtain a twofold cycle as op- 
posed to a fourfold one. Hence, the number of echoes in a 
cycle is observer-dependent but the existence of a cyclic 
feature is universal and associated with the trapping of 
energy near the black hole. 

Consider the top panel of Fig. [13] The first pulse is 
formed by the direct signal and the N = 1 caustic echo. 
The next pulse arrives after a full revolution time Tf u n 
and consists of the N = 2, 3 caustic echoes. (Compare 
Fig.[9]or Figs. [3][6] with an observer along the caustic line 
to visualize the simultaneous arrival of echoes.) 
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FIG. 12: Top: The caustic echoes of a square pulse source 
lasting for 6A/ and located at 6M on the positive x axis as 
seen by an observer on the positive z axis at future null infin- 
ity. Bottom: The direct signal of the square pulse (black), its 
Hilbert transform (blue), and the first caustic echo appropri- 
ately shifted and rescaled to match the amplitude of the blue 
curve. Despite the long time scale of the square pulse, its cor- 
responding wavefronts are still Hilbert transformed through 
caustics. 



The echoes on the caustic line also undergo a Hilbert 
transformation. The middle panel of Fig. [13] shows 
the first (N = 0,1) and the second (N = 2,3) de- 
generate echoes for an observer on the negative x axis 
(7 = 7r), while the bottom panel of Fig. Il3] shows the 
first (N — 1,2) and the second (N — 3, 4) degenerate 
echoes for an observer on the positive x axis (7 = 0). 
The N = 0, 1 degeneracy for the 7 = it observer leads 
with a Gaussian profile which is Hilbert transformed at 
its peak to a Dawsonian for the trailing half of the signal. 
The N — 1,2 degeneracy for the 7 = observer leads 
with a Dawsonian profile which is Hilbert transformed at 
its zero crossing to a Gaussian for the trailing half of the 
signal. The next degenerate echoes are the Hilbert trans- 
forms of the previous ones as indicated in the middle and 
bottom right panels of Fig. [l3j The red curves are the 
Hilbert transforms of the previous degenerate echoes to 
the left, suitably shifted in time and rescaled in ampli- 
tude for comparison. 

The amplitude of the wave becomes singular along 
the caustic line in the geometrical optics approxima- 
tion. However, given that the echoes can be charac- 
terized via Hilbert transforms whenever the wavefront 
passes through a caustic, the time dependence of the sig- 
nal may be captured by the geometrical optics approxi- 
mation. To test this, we fit the data to the profile implied 
by the convolution of the Green function in the geomet- 
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FIG. 13: Top: The evolution of the scalar field as measured 
by an observer at future null infinity on the negative x axis. 
Middle: Degenerate echoes as measured by the same ob- 
server. Because the source is also on the x axis, the observer 
only sees a twofold cycle. Also shown is the Hilbert transform 
of the signal (red) in the left panel (rescaled and shifted to 
compare with the caustic echo) and the free fit to the geo- 
metrical optics approximation (dashed). Bottom: Same as 
above except that the observer is on the positive x axis. Note 
the differences in the caustic echo profiles. 



rical optics approximation ^ with the source (J2]>- More 
specifically, the N = 0, 1 degeneracy for the 7 = ir ob- 
server is fit to 



!e -(t-6)V(2d 2 ) +c 0( t -b)D{t-b,d) 



(11) 



where b is approximately the time of arrival To = X\, 
while for the N = 1, 2 degeneracy for the 7 = observer 
the fit is 



aD(t-b,d) + c9(t-b) e -(*^) 2 /(2d 2 ) 



(12) 



with b approximately T\ = T2. In both cases, the am- 
plitudes (a, c), the time of arrival (6), and the width are 
arbitrary (d). The results of this free fit are shown in 



the middle and bottom right panels of Fig. 13 (dashed 
lines). Despite its breakdown, the geometrical optics ap- 
proximation captures the basic time-dependence of the 
profiles fairly well. 

The wavefront is focused and its amplitude is magnified 
along the caustic line (see Figs. 3|6). We next character- 
ize this magnification at the caustic relative to the echoes 
it generates in relation to the scale of the wavefront. The 
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computations of such relations go beyond the geometrical 
optics approximation and require elements of geometric 
diffraction theory. In [53] Kay and Keller show that the 
peak amplitude at a caustic with a perfect focus increases 
as (J -1 / 2 , where a provides the scale for the width of the 
wavefront. 

Note that the caustic and its echoes have different 
shapes. Therefore, a comparison of their peak amplitude 
might be misleading. Instead, we compare the energy 
radiated to an observer at infinity during a certain time 
interval /, 

E=\Jdt\4>{t)\\ (13) 

To make sure that the amplification factor is not in- 
fluenced by the shape of the echoes, we compute it with 
respect to both Dawsonian and Gaussian shapes. The 
compared signals are plotted in the top panel of Fig. [14] 
We label the measurement of the caustic (N = 1,2) by 
an observer at 7 = as A. An observer at 7 = it/2 
measures not the caustic but its echoes. The Dawsonian 
echo (N = 1) is labeled as B, the Gaussian one (N = 2) 
as C. 



We plot the amplification on a log-log scale as a func- 
tion of Gaussian widths a in the bottom two panels of 
Fig. [14] The amplification factor can be fit accurately by 
an inverse power law of the form aa p where (o.b,Pb) = 
(23.7,-1.03) (bottom left) and (a c ,pc) = (26.1,-1.02) 
(bottom right). Hence, the data suggests that the energy 
amplification at the caustic goes as o~~ x as the width of 
the Gaussian source becomes more narrow. This seems 
in accordance with the calculation of Kay and Keller for 
caustics with a perfect focus [23] , 

We emphasize, however, that we cover only a small 
range of scales [a E [0.1,0.6]) and that we did not math- 
ematically derive the scale dependence of the energy at 
caustics formed by a Schwarzschild black hole. Includ- 
ing smaller values for a will most likely follow an inverse 
power law but may change the value in the fit. Therefore, 
our numerical result should be interpreted carefully. 

There are extensive studies on the detailed behavior of 
solutions to wave equations near caustics, including not 
only their energy magnification but also their shapes in 
different directions away from the caustic line. Such an 
analysis is outside the scope of this paper but provides 
an interesting subject for future research. 
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FIG. 14: Top: A caustic (A) at 7 = (black), the N = 1 echo 
(B), and the N = 2 (C) echo along 7 = it/2 (red). Bottom: 
The amplification of the energy of A relative to B (left) and 
C (right) along with an inverse power law fit (dashed). 

We compute the energy integral over a sufficiently large 
time interval centered at the pulse so that the energy 
difference to a larger interval is negligible (less than 1%). 
The amplification is then computed via 



Amplification 



E A 
Eb.c 



(14) 



E. Propagation within the light cone 

We have shown that essential properties of the Green 
function can already be understood in the geometrical 
optics limit. We also discussed effects due to the finite 
wavelength of the wave, which plays a role in the energy 
magnification at caustics. In this section, we consider 
features at the largest scale, the spacetime curvature. 

It is well-known that the Green function propagates 
also inside the null cone in curved spacetimes due to 
the backscatter of the field off the background curvature. 
The field within the null cone decays polynomially fol- 
lowing the so-called Price power law [24]. The rate of 
the decay for idealized observers at null infinity is dif- 
ferent than the rate at finite distances [53]. While the 
null infinity rate is the relevant one with respect to ide- 
alized observers, the finite distance rate is the one that 
effects the self-force on a particle. Both rates are plotted 
in Fig. [15] 

The theoretical decay rates are valid asymptotically in 
time. One can obtain a good estimate of the rates by 
computing a local in time rate along the worldline of an 
observer at r — R defined by 



PR(t) 



d]n\<f>(t, R)\ 
dint ' 



(15) 



The theoretical asymptotic decay rate for a generic scalar 
perturbation in Schwarzschild spacetime is for idealized 
observers at null infinity p^ = — 2 and for finite distance 
observers pn = —3. Figure [15] confirms the approach of 
the local power to the asymptotic rate. The theoretical 
predictions are plotted as dashed lines. 
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FIG. 15: Polynomial decay rates asymptotically approaching 
p — —2 for an idealized observer at null infinity (top solid 
curve), and p = —3 for a finite distance observer at r = 4M 
(bottom solid curve). The theoretical asymptotic decay rates 
are plotted as dashed lines. Caustic echoes are still visible in 
the local power plot at late times. 
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FIG. 16: The visibility of a caustic echo in the local power 
at null infinity for three different values of the width of the 
initial Gaussian source Q. More caustic echoes are visible at 
late times for narrower Gaussians. 



In discussions of gravitational lensing it is often stated 
that, under certain assumptions and with ideal instru- 
ments, one would see infinitely many images of a source 
with an eternal worldline |26) . Considering the geomet- 
rical optics limit, one might be inclined to think that for 
arbitrarily narrow Gaussians we would observe arbitrar- 
ily many caustic echoes. We show in Fig. [16] a caustic 
echo visible in the local power of the signal at null infin- 
ity for three different values of a for the Gaussian source 
([2]). The caustic echo is more pronounced for smaller 
values of a indicating that the narrower our source, the 
more caustic echoes we see in the late time signal. 

However, it is also clear that for any positive value of 
a we can only see a finite number of caustic echoes be- 
cause the amplitude of each echo decays exponentially in 
time whereas the backscatter off curvature decays only 
polynomially. This observation implies that, even with 
ideal instruments, we can measure only a finite number 
of "images" for any given source because eventually the 
polynomial decay will win over the exponentially weak- 



ening echo. Proofs of infinitely many images of sources 
typically consider the geodesic structure of spacetimes 
and neglect backscatter [27] . 



F. Putting the pieces together 

We presented in previous sections a quantitative un- 
derstanding of each feature in the evolution of a scalar 
perturbation triggered by a narrow Gaussian wave pack- 
age. Dissecting the signal observed at infinity (Fig. [2]) 
we explained the arrival times, the exponential decay, 
the shapes of the echoes, and the late time behavior of 
the scalar field. In this section, we bring these elements 
together into a heuristic formula that captures the essen- 
tial features of the evolution remarkably well. 

The shape, the arrival times, and the exponential decay 
of caustic echoes can be accurately described in the geo- 
metrical optics limit as discussed in Sees. |IIIB| and [TTl C| 
We derived an analytic expression for the Schwarzschild 
Green function in this limit (Appendix |B|) , presented for 
the first time in this paper as far as we are aware. Using 
this analytic knowledge we approximate the signal by 



^geom — Ao 



e -(t-T a f/(2al) +e -\{t-T ) 



N 



irvK 1 2 



n=l 



X e -(t-T n -STf/(2al) _ w{t _ ^ _ ^ 



(16) 

where A 0.096225 is the Lyapunov exponent for unsta- 
ble photon orbits, To is the time of arrival of the direct 
signal to the observer, and T„ is the time of arrival for 
the n th caustic echo (with T_i = 0). The latter is com- 
puted from T n _i + 5T + AT n (-f) using The quantity 
ST accounts for the fact that the time between succes- 
sive caustic echoes is not immediately given by AT n (-j) 
(Fig. [7]). We found it sufficient to approximate ST by 
T 2 — Tq — Tf u ii, which is independent of 7. Only T and 
ST (or T2) are extracted from the numerical simulation. 
We set Aq — 0.016, a a = 0.28, and the number of caustic 
echoes as N — 15 specifically for the observer of Fig. [2j 

The geometrical optics formula for the scalar field 
agrees with the numerical data in the early part of the 
signal (top panel of Fig. 17 1. To match the late time be- 
havior of the field we use the polynomial decay discussed 
in Sec. |III E| We include the tail only after the direct 
signal has reached the observer and set 



<kaa = Cf0(t-r o ), 



(17) 



with C = —0.06 and p = —2.25. Note that the asymp- 
totic decay rate would be p = —2. The difference is due 
to the high order correction terms which are neglected in 
the above formula. One can include such correction terms 
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FIG. 17: Reconstructing the numerical solution (black) using 
heuristic formulae (red). The top panel shows the structure 
of the field from the geometrical optics limit. The bottom 
panel includes the power law tail contribution. Despite the 
simple nature of the approximations, the field's rich dynamics 
are captured remarkably well. 



by constructing more general fitting functions (see for ex- 
ample [Mj)- Recently, Casals and Ottewill computed the 
first three orders in the power law for the Schwarzschild 
Green function [55], which can also be used to improve 
the description of the tail. We found the simple expres- 
sion above sufficiently accurate for our purposes. Better 
expressions for the tail fit might be needed for compar- 
isons with more accurate approximations to the Green 
function, or for self-force computations. 

Our heuristic formula for the full field is the sum of the 



geometrical optics contribution ( 16 1 and the tail contri- 
bution (17 1. Considering the simplicity of the assump- 
tions that go into our heuristic expression, its agree- 
ment with numerical data is remarkable (bottom panel 
of Fig. ~ 



17) 



Note that the addition of the tail not only improves 
the agreement with the data at late times but also cap- 
tures the zero crossings between the caustic echoes at 
t/M s» 57 and 86. In typical numerical studies of the 
tail, the polynomially decaying part of the signal plays 
a role only at late times. Here, we see that it improves 
the fit between the echoes at relatively early times. This 
is because backscattering contributes even right after the 
direct signal has passed the observer. This feature is 
well-known and is captured in Hadamard's ansatz for the 
Green function whenever x and x' are connected by a 
unique geodesic [301 131] • 

Our heuristic formula should fail along the caustic line 
(the x axis in our simulations) because the geometrical 
optics approximation breaks down at caustics. We find 
that ( |l6[ ) and (17) on the caustic line works well for the 
high-frequency part of the signal (i.e., the caustic echoes 



themselves) and for describing the late-time tail (Fig. 18 ), 
but cannot capture the effects of early-time backscatter- 
ing (compare to Fig. 17 for the observer at 7 = 7r/2). 




200 



FIG. 18: Reconstructing the numerical solution (black) as 
seen by an observer at null infinity on the positive x axis along 
the caustic line. The analytical approximation (red) disagrees 
significantly with the data between the caustic echoes at early 
times. 



cur at the correct times. Near a caustic the wavefronts 
are focused and magnified thus distorting the shape of 
the wavepacket (see Figs. 3 6). Hence, the echoes have 
a comparatively larger amplitude than the direct signal 
(Fig. 18). These features are not accounted for in the 
geometrical optics limit and require a more complete de- 
scription of wave optics in curved spacetime. 

Further analysis is needed to improve on our heuristic 
formula. For example, one can include corrections to the 
geometrical optics approximation using the expansions 
presented in [8] or combine other quasilocal expansions 
of the Green function [32] and its Pade resummation [J] . 
These potential directions for research are left for future 
work. 



Also, the formula has zero crossings which do not oc- 



IV. CONCLUSIONS 

We numerically solved the scalar wave equation on 
a Schwarzschild background with a narrow Gaus- 
sian source |2]) thus yielding an approximation to the 
Schwarzschild Green function ([3]). For the numerical 
computations we used the Spectral Einstein Code SpEC 

The main result of our study is the numerical approx- 
imation of the full retarded Schwarzschild Green func- 
tion in the time domain. This computation allowed us 
to demonstrate a cyclic singularity structure due to the 
trapping of energy near the photon sphere. 

The numerical evolution, visualized in Figs. 3][6 and 
in [13], proceeds as follows. The narrow Gaussian 
source triggers a high-frequency wavefront that propa- 
gates along null geodesies in accordance with geometri- 
cal optics. Part of the energy of the initial wavefront 
gets trapped at the black hole horizon and leaks out to 
infinity decaying exponentially in time. At each half rev- 
olution around the photon sphere the trapped wavefront 
forms a caustic resulting in echoes that propagate out 
to infinity. The wavefront undergoes a Hilbert transform 
at the caustics causing a —tt/2 phase shift in its profile, 
which results in a fourfold cycle of caustic echoes seen 
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by generic observers. This cycle is known as the four- 
fold singularity structure of the retarded Green function 

[DHHZIISI. 

The fourfold structure is a recent discovery in the con- 
text of curved spacetimes but, in hindsight, it is a simple 
consequence of well-known facts. The 7r/2 phase shift 
of wavefronts at caustics has already been discovered in 
the late 19th century and is widely known as the Gouy 
phase shift in optics [20] [3T] . Combined with the notion 
of trapping (e.g., due to the presence of a photon sphere), 
this knowledge immediately implies a fourfold cycle in a 
spherically symmetric, asymptotically flat spacetime 1 . 
The viewpoint of a phase shift of the Green function at 
caustics seems more fundamental than the fourfold struc- 
ture. For example, a phase shift may be observed in re- 
gions of moderate curvature due to gravitational lensing 
where no trapping occurs, whereas trapping is required 
for the observation of n-fold cycles. In addition, a non- 
generic set of observers, those on the caustic line, see a 
twofold cycle. Hence, the number of echoes in a cycle is 
observer-dependent. 

With the full numerical approximation to the Green 
function at hand, we also obtained an interesting result 
regarding gravitational lensing. If backscatter off back- 
ground curvature is included, only a finite number of 
images are visible, even with ideal instruments, for any 
source of finite wavelength because caustic echoes decay 
exponentially in amplitude while backscatter (typically 
neglected in lensing studies) decays only polynomially. 

Further results can be summarized as follows. The 
arrival and decay of successive caustic echoes are con- 
sistent with the orbital period and Lyapunov exponent 
of null geodesies trapped at the photon sphere. The ar- 
rival times of successive half-period echoes depend on the 
source-observer configuration, with degeneracy resulting 
in a twofold cycle when the source and the observer are 
aligned with the black hole. This is the caustic line along 
which the amplified signal propagates. The energy mag- 
nification at the caustic follows an inverse power law with 
the scale of the wavefront in accordance with predic- 
tions going beyond the geometrical optics limit [53] . The 
backscatter follows well-known polynomial decay rates 
at null infinity and at finite distances, which limits the 
number of echoes measurable by observers. 

We formulated an analytic approximation to the 
Schwarzschild Green function in the geometrical optics 
limit including wavefront propagation through caustics. 
Combined with our quantitative understanding of the dy- 
namics, this formula allowed us to capture the essential 
features of the observed signal in a heuristic expression. 
Given the simplicity of our approximations, it is remark- 
able that the expression agrees so well with the numerical 



calculation. 

The only sources of error in the numerical computa- 
tion of the Schwarzschild Green function are the finite 
width of the Gaussian source ^ and the truncation er- 
ror due to discretization of the scalar wave equation. The 
boundary error typically present in such simulations is 
eliminated by hyperboloidal scri-fixing in a layer |lll I12j . 
We emphasize that hyperboloidal compactification is not 
necessary for the numerical study of the Green func- 
tion near the black hole. Many of our results are re- 
producible using standard foliations. However, hyper- 
boloidal compactification helps focusing the resolution 
to the strong field domain without contamination from 
artificial boundaries and without wasting computational 
resources on the asymptotic domain. It allows us to com- 
pute the long-time evolution accurately at low cost, and 
provides us direct numerical access to measurements of 
an idealized observer at future null infinity. 

The error scales in our computation are related. To 
evolve a narrow Gaussian source, we need a large num- 
ber of collocation points for the spectral expansion of the 
variables. Failure to do so results in Gibbs phenomena 
and contaminates the evolution. (Some high frequency 
noise in the late-time evolution can be seen in the lower 
curve of Fig. 15). Numerical convergence tests and the 



1 Trapping may also occur globally in cosmological spacetimes 
(e.g., the Einstein static universe), which may result in a dif- 
ferent phase shift and therefore a different n-fold structure. 



agreement between theory and experiment indicate that 
our errors are small but we did not present a detailed 
error analysis in this paper. Such an analysis will be 
required when comparing the numerical simulations to 
accurate approximations of the Green function or when 
computing the self-force acting on a particle using the 
numerically computed Green function. For these calcu- 
lations, horizon-source ratios larger than 10:1 might be 
necessary. The numerical method must be further im- 
proved for very high ratios, possibly using a second or- 
der in space formulation [33], implicit-explicit time step- 
ping [34], and adaptive mesh refinement. There are also 
adapted methods to solve for high-frequency wave prop- 
agation, such as the frozen Gaussian beam method |35j 
or the butterfly algorithm [36], that may improve the 
efficiency of the numerical simulation considerably. 

An immediate extension of our study would be the nu- 
merical computation of retarded Green functions in Kerr 
spacetime. In rotating black hole spacetimes, frame drag- 
ging plays a role in the wavefront propagation. There is 
no photon sphere in Kerr spacetime; instead, there are 
spherical photon orbits bounded by the location of pro- 
and retrograde circular photon orbits [37, which may par- 
ticipate in the Cauchy evolution in a similar way as does 
the photon sphere in Schwarzschild spacetime. The large 
I limit of the quasinormal mode spectrum and its rela- 
tion to spherical photon orbits in Kerr spacetime have 
recently been analyzed in [35]. Such analytic knowledge 
can be combined with numerical experiments to reveal 
the structure of the Green function in Kerr spacetimes 
and to provide good approximations to it. The visualiza- 
tion of a numerical Kerr Green function simulation using 
SpEC can be found in |3"9"] . 
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An interesting application of numerical approxima- 
tions of Green functions would be the computation of 
the self-force on a small compact object moving in a su- 
permassive black hole spacetime. Note that the Green 
function provides a very general way of solving for wave 
propagation in a background spacetime. If good approx- 
imations to it can be found, possibly augmented by data 
from numerical computations, the method can be applied 
to essentially any problem on the given background via 
suitable convolutions. To this end, it may be useful to 
have sparse representations of the numerical Green func- 
tion. 

The improvement of Green function approximations 
may benefit from developments in various fields. Forma- 
tion of caustics in black hole spacetimes is an example 
of the emergence of discrete structures from continuous 
ones, a hallmark of catastrophe theory [40]. A develop- 
ment of geometric diffraction theory and wave optics for 
black holes may allow us to compute the Green function 
through caustics without matching arguments while also 
providing more accurate analytical approximations. Ex- 
citing developments in these directions can be expected 
in the near future. 
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where da 2 = ef$ 2 +sin d 2 dip 2 is the standard metric on the 
unit sphere. Ingoing Eddington-Finkelstein coordinates 
are constructed from Schwarzschild coordinates by the 
time transformation t = ts + 2M\n(r — 2M), which leads 
to the metric 
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+ (l + 2 M\ d r 2 +r 2 da 2 . (Al) 



The metric is regular at the horizon, r = 2M, because 
the slices of the iEF foliation penetrate the future horizon 
without intersecting each other, as opposed to the slices 
of the standard Schwarzschild foliation which intersect 
at the bifurcation sphere. In a similar construction, iEF 
slices that intersect at spatial infinity can be transformed 
to hyperboloidal slices that foliate null infinity [4Tj . The 
hyperboloidal layer method that leads to such a slicing 
consists in the introduction of a compactifying coordinate 
in combination with a suitable time transformation. The 
compactifying coordinate p can be defined via 
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Appendix A: Hyperboloidal layer with excision 

In our numerical simulations we use hyperboloidal scri- 
fixing [TT] in a layer [T^] attached to a finite Schwarzschild 
domain in ingoing Eddington-Finkelstein (iEF) coordi- 
nates. In this Appendix we describe the construction of 
such coordinates. 

The Schwarzschild metric in standard Schwarzschild 
coordinates {ts, r, ip} reads 

* = -(l--J*S+(l--J dr 2 + r 2 d* 2 , 



r=£, where Q = 1 - (j^j) ~ Z ) • 

(A2) 

Here, is the Heaviside step function. The compact- 
ifying layer is attached at the interface p — I and has 
a thickness of S — I in local coordinates. The zero set 
of f2, namely p = S, corresponds to infinity. To avoid 
loss of resolution near infinity, we combine this spatial 
compactification with a time transformation. The time 
transformation is constructed requiring invariance of out- 
going characteristic speeds [HI 22], or equivalency, of the 
outgoing null surfaces in local coordinates [121 21] . Out- 
going null surfaces satisfy u = t — (r + AM \og(r — 2M)) . 
The invariance condition translates to 



t - (r + 4M log(r - 2M)) = r-{p + AM log(p - 2M)) . 



This relation defines the hyperboloidal time coordinate 
r used in the compactifying layer. The metric resulting 
from these transformations is singular up to a rescaling. 
The conformally rescaled metric reads 
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where L = Q — p dQ /dp. The metric looks rather com- 
plicated but notice that for Q = 1 (or equivalently, for 
p < I) it reduces to the iEF metric (Al) and at null 
infinity, Q = 0, it takes a purely outgoing form. As a 
consequence, characteristic speeds incoming with respect 
to the bulk vanish at the boundaries. 



singularity 



a 0.0 





FIG. 19: Radial characteristic speeds in ingoing Eddington- 
Finkelstein coordinates with a hyperboloidal layer attached. 
The dashed vertical red line depicts the interface to the layer 
at 7 = 8. Null infinity is at the coordinate location S = 12. 
These parameters are also used in the numerical simulations. 



FIG. 20: Conformal diagram of an ingoing Eddington- 
Finkelstein foliation with a hyperboloidal layer. The dashed 
red line corresponds to the interface. For visualization we 
have chosen 7 = 3 and S = 4. For the numerical computa- 
tions we set 7 = 8 and S = 12 as in Fig. |19| 



The radial characteristic speeds read 

_ p~2M 
C+ ~ p + 2M' 

2 (p- 2M)(p- 2MQ) 

2Lp{p - 2M) + n 2 (p + 2M)(p - 2Mf2) ' 

The outgoing speed c+ has the same form as in iEF co- 
ordinates. Consequently, it vanishes at the event horizon 
r = p = 2M. Similarly, the ingoing speed c_ vanishes at 
null infinity, p = S where fl = 0. Figure [19] shows the ra- 
dial characteristic speeds in these coordinates. The con- 
formal diagram of the iEF foliation with a hyperboloidal 
layer is depicted in Fig. [20} 

In the hyperboloidal layer we solve the conformally 
transformed scalar wave equation ([!]) that reads [IS] 

□ -R$ = Sl- 3 S(x), (A3) 
6 

where $ = 0/0 and R is the Ricci scalar of the conformal 
background. The divisions by the conformal factor do 
not pose any difficulty because the scalar field falls off 
suitably and the source has compact support. 



Appendix B: Green function in geometrical optics 

The numerical solutions presented in Sec. |III| can be 
regarded as approximations to the Green function for 
the scalar wave equation (nj). We derive here, for the 
first time to our knowledge, the scalar Green function in 
Schwarzschild spacetime in the geometrical optics limit, 
which is used in the main text to compare numerical data 
with analytical approximations and to explain the four- 
fold cycle and the profiles of the caustic echoes. See [5T] 
for a review of the geometrical optics limit in the field of 
optics. 

In Schwarzschild spacetime there exists a time-like 
Killing vector, d t , that allows for the Green function to 
be decomposed into frequency modes via 

/OO 1 
^e-^G^xy). (Bl) 

In the frequency domain, using coordinates in which 
g u = 0, the Green function equation ^ becomes a 
Hclmholtz equation 

AG(w; x, x') + w V (x)G(w; x, x') = 0, (B2) 

where A = g^VjVj is the curved space Laplacian (with 
i,j = 1,2,3), and we assume that x ^ x' so that the 
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FIG. 21: The retarded Green function as a mapping from a 
wavefront with phase ip' to a wavefront with phase ip > ip' . 



Dirac delta distribution in ^ vanishes. The function 
it(x) is defined by 



u 2 (x) = - 5 «(x). 



(B3) 



2M 



In standard Schwarzschild coordinates m 2 (x) = 1 

In the geometrical optics limit, the frequency of the 
mode approaches infinity (the wavelength of the mode is 
much smaller than the background curvature scale). We 
make the following ansatz for the asymptotic expansion 
of the Green function, 



G( W ;x,x') 



x') eiwT ( XiX ' 

n=0 



(iu)) n 



(B4) 



Here, T(x,x') is called the eikonal and its level surfaces 
are surfaces of constant phase, called wavefronts. The 
covector V^T is everywhere normal to this surface. Note 
that this ansatz is not covariant as time is a preferred 
direction. This may seem undesirable but our static ob- 
servers measuring the field, as in Fig. [2j are Killing ob- 
servers. Additionally, the geometrical optics approxima- 
tion emphasizes time over space via the assumption of 
only high-frequencies contributing to the field. Hence, 
the lack of manifest covariance is suitable for our pur- 
poses. 

Different level surfaces correspond to wavefronts with 
different phases of the Green function. We may thus 
parameterize these wavefronts by the phase, which is an 
affine parameter denoted by (p. As the field increases 
phase, and thus evolves forward in time, one moves from 
one wavefront to the next. This motion can be regarded 
as a mapping that allows one to identify paths, or rays, 
that connect specific points on different wavefronts. The 
Green function in ( B4 ) thus causally propagates points 



from one wavefront to another (Fig. 21 ). This viewpoint 
will be useful below. 



Substituting the leading order term of the sum in ( B4 ) 



into (B2) and equating like powers of u> gives the fol- 



lowing set of equations for the eikonal T(x, x') and the 
leading order amplitude j4o( x j x ')j 



V,TV ! T = u 2 (x) 
2ViA V*T = -A Q AT. 



(B5) 
(B6) 



The first equation is called the eikonal equation and the 
second is called the transport equation. We will solve 
these equations in turn. 



The solution to the eikonal equation ( B5 ) can be found 
by first defining a momentum vector pi = VjT(x), which 
is normal to the wavefront so that the eikonal equation 
reads p 2 = it 2 (x). This defines a Hamiltonian for the 
rays connecting one wavefront to another 



H 



~(p iPj g ij (x)-u 2 {x)). 



(B7) 



The Hamiltonian vanishes when the eikonal equation is 
satisfied. Hamilton's equations give 



p l = V i T(x) 



dx 1 
dip 



1. 



(B8) 
(B9) 



where D/dip — (dx 1 / dip)V i is the covariant parameter 
derivative along a trajectory with coordinates x 1 (ip) (i.e., 
a ray). Solving these ray equations thus amounts to solv- 
ing the eikonal equation. 

The square of the ray's velocity is 



dx l dxi 
dip dip 



ViTV l T = m 2 (x) = -<?" (x). (BIO) 



Collecting all terms on one side and multiplying by dip 2 
yields 



g t1: (x)dip 2 + g ij (x)dx t dxj = 0. 



(Bll) 



Identifying dx° = dt = u 2 (x.)dip, we see that the rays in 
the geometrical optics limit are null because dx^dx^ = 
0. The solution to the eikonal equation thus amounts 
to solving for the trajectories of the null rays, which is 
equivalent to solving the geodesic equation for two light- 
like separated points, one on each wavefront. 

Next, we solve the transport equation (B6). Divid- 
ing both sides by Aq, noti ng that AT = ViX l (ip), and 
V*TV 2 = D/dip (from jB8|) yields 



din A Q 

dip 



(B12) 



If an infinitesimal area on the wavefront centered at 
x' = x(ip') is evolved along the null geodesies to phase 
tp > ip' then the subsequent evolution can be viewed as a 
change of variable since x = x(y>) depends on the "initial 
data" at x' and implies the coordinate transformation 
x = x(tp;x!). The Jacobian relating these areas is 



J(x,x') 



d(x° ■ ■ ■ x 3 ) 
d(x'°---x' 3 ) 



detM 



Tr In M 



(B13) 



where the matrix M has components M l j = dx l jdx'K 
Taking the parameter derivative of the Jacobian gives 



dJ _ jdx H d±i 
dip dxi dx n 



= J 



.dx 11 dxi g x k 
dx^ Q x k dx' 



= J 



.dx 1 
dx 1 



(B14) 
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Using diX 1 — Vii J — T\ a x a and 



dip 



(B15) 



where h is the determinant of the spatial metric g^, we 



obtain from (B12 1 



d In A n 2 „ . ,■ d / T , 

dip dip 



(B16) 



This has the solution 
Ao(x,x') = A (x',x") 



J(x',x") y/MxQ 

j( x ,x') y^xj 



1/2 



, (B17) 



provided that x = x((p) has not passed through any caus- 
tics and where j4o(x',x") is initial amplitude (x" is ei- 
ther a point on an earlier wavefront or can be equal to 
x' in which case J(x',x') = 1). Hence, the amplitude 
of the wave at phase ip 1 for a given infinitesimal area 
around the point x' = x(y>') on the wavefront will evolve 
along the null ray to x = x(<p) with phase ip > ip' ac- 
cording to (B17). Notice that the amplitude at x de- 



pends on the expansion or contraction of the infinitesi- 
mal area element on the wavefront during its evolution. 
This is equivalently described by the expansion or con- 
traction of a congruence of null geodesies connecting x' 
and x. When the ray passes through a caustic the Jaco- 
bian J(x, x') changes sign relative to J(x',x"). We will 
say more about this in Appendix [C] 

Assuming that the ray has not passed through any 
caustics, we may reconstruct the time-domain Green 
function in (Bl) using (B4) and ( B17 1 to find at lead- 
ing order in 1/lu, 



G{x,x') 



A,(x,x') 
x Re 



J(x', X ") y/^xQ 



1/2 



du 
2^ 



-iw(t-t'-T(x,x')) 



(B18) 



where we take the real part of the integral because the 
field in the time domain is real valued. The evaluation 
of the integral gives 



G(x, x') oc S(t — t' — T(x, x')j . 



(B19) 



Hence, the geometrical optics approximation yields a con- 
tribution to the Green function whenever the time differ- 
ence between x' and x satisfies t — t' = T(x, x'), which is 
just the time it takes for a null geodesic to connect these 
two points. 

Finally, we note that the remaining amplitudes 
yl„(x, x') in ( |B4[ ) can be determined through a recursion 
relation (for n > 0), 



2V l TV^„+i + A n+1 AT = -AA r . 



(B20) 



which follows from substituting (B4) into (B2) and set- 
ting the coeffici ents o f each power of lj to zero. With 
A known from (B17) and the rays found by solving the 



null geodesic equation, we have the Schwarzschild Green 
function in the geometrical optics limit. The rays are 
null at all orders of 1/u). Hence, the geometrical optics 
limit is inapplicable to the propagation into the null cone 
of x' and thus ignores backscattering off the background 
curvature. 



Appendix C: Hilbert transform through a caustic 



In Appendix |B| we solved the transport equation ( B6 ) 
for the field's amplitude in the geometrical optics limit. 
The solution (B17) assumes that the quantity under the 
radical is positive. However, if x = x(ip) evolves through 
a caustic then the Jacobian (B13 ) changes sign. We show 



in Fig. [22] a cartoon of a wavefront area element as it 
passes through a generic caustic. 





Caustic 




7' 













FIG. 22: A cartoon of the change in orientation of an in- 
finitesimal area element on the wavefront's spatial section as 
it passes through a caustic in Schwarzschild spacetime. Only 
one dimension of the area element is compressed. 



The sign change of the Jacobian throu gh a c austic im- 
plies a purely imaginary amplitude in (B17). Hence, 



the amplitude changes phase by ±7r/2 through a caus- 
tic (20] [21] . The sign of the phase shift is determined by 
matching to the full solution from numerical simulations. 
It cannot be computed within the geometrical optics ap- 
proximation because the latter breaks down at caustics 
(where J(x,x') = and A(x, x') becomes singular). 

Aft er some trial and error, we find that the phase of 
(B17) when passing through one caustic is given by 



exp 



-sgn(w) 



(CI) 



For positive frequencies the phase is — i while for nega- 
tive frequencies it is +i. The frequency-dependence in 
(CI ) is needed for the reconstruction of the time-domain 



Green function in the geometrical optics limit. Comput- 



ing (Bl ), upon substituting the leading order term of the 
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1 /oj expansion from ( B4 ) , gives 

J(x',x") 



G(x,x')~\A (x', X ") 



J(x,x') 
x e(t-t' -T (x,x')) 
x Re 



I hp) 



1/2 



dio 
2^ 



g-iw(t-i'-Ti(x,x')) e -^sgn(w) 



(C2) 



where we have explicitly taken the real part as in 
T/v(x,x'), with A" = 0,1, is the eikonal for the null 
geodesic connecting x and x' and passing through AT 
caustics, and the step function ensures causality (namely, 
that the T\ eikonal is larger than To)- The frequency in- 
tegral gives 



(C3) 



(t-f -TiCx.x')) 



As discussed in Sec. |III C| when the Green function is 
convolved with the source S(x) from (|| one finds that 
passing through one caustic transforms the field from a 
Gaussian profile to that of a Dawson's integral (or "Daw- 
sonian" ) , which is the convolution of ( C3 1 with the source 
S(x) [32] . Hence, the choi ce gi ven in ( |C1[ ) for the sign 
of the root of j4q(x,x') in (B17) is justified through this 
matching calculation. 



Comparing (C3| to (B19) we note that the latter has 



support only along the null ray connecting x and x' 
while the former has a much broader domain of sup- 
port. Both expressions are singular on null rays (where 
t-t' = T (x,x') or r!(x,x')), but (pjj) does not vanish 



in the forward light cone thus adding to the tail. Passing 
through one caustic causes the sharp delta distribution 
to become smeared out into the future null cone of x' 
beyond the caustic. 



The frequency integral in ( C2 I is the Hilbert transform 
of 8(t - t' - Ti(x,x')). For a function /(i), the Hilbert 
transform in the frequency domain is defined as 



H[f(t)} = 



duj 
2n 



/Me" 



T/2sgn(w) 



(C4) 



where /(w) is the Fourier transform of f(t). The effect 
of the Hilbert transform is to shift the phase of f(tu) by 
— 7r/2. Two applications of the Hilbert transform gives 
H 2 [f(t)} = H[H[f(t)}] = -/(*), which is a -tt phase 
shift of f(t). Four applications of the Hilbert transform 
gives back the function itself, 



(C5) 



Thus, the Hilbert transform explains the fourfold cycle 
of the retarded Green function through caustics. This 
allows us to generalize ( |C2[ ) to include all the null rays 
that connect x and x' (of which there are infinitely many 



V- G(w;£l,b') -A |- G(u;x,x 2 ) -| 

|— G(u;xi,x') — |— G(u>;x2,£i) — j 






FIG. 23: Top: The retarded Green function in the geomet- 
rical optics limit for a null ray connecting x' to x through a 
single caustic. The geometrical optics approximation breaks 
down at the caustic (yellow star) but is valid from the wave- 
fronts with phase <p' to tpi > ip' and from the wavefront with 
phase <p2 to tp > <j>2 . Bottom: In Maslov's theory the caustic 
is displaced by a phase space canonical transformation yield- 
ing a related Green function that allows for propagating the 
wavefront with phase tp' directly to tp2 via <pi . 



for a Schwarzschild black hole [2"T]h 

J(x',x") 



G(x,a/)~ |A)(x',x" 



N=0 




J(x,x') 

/OO J 

Sgn(u) i ( C6 ) 



x e - iw (t-t'-T,v(x,x')) e - i ^ 



where T_i(x, x') = 0. This formula can be written more 
succinctly as 

oo 

G(x,x')~ J2 lA^xOl^-t'-lV-i^x')) 



N=Q 



Re{H N [5(t-t' -T N (x,x'))]}, 



(C7) 

where Tjv(x, x') is the eikonal for the null geodesic that 
passes through A^ caustics and Hm denotes the Hilbert 
transform applied N times to its argument. 

As mentioned at the beginning of this Appendix, the 
geometrical optics approximation breaks down at caus- 
tics where many points of a wavefront get mapped to 
the same (caustic) point. The Jacobian J(x, x') van- 
ishes and the amplitude of the field Ao(x, x') diverges 
as J(x,x') -1 / 2 . The geometrical optics approximation is 
valid only between successive caustics. One may use the 
geometrical optics Green function to evolve wavefronts 
before and after a caustic but not through (Fig. 23). 

Remarkably, the breakdown of geometrical optics at 
a caustic can be averted (or rather, displaced) using 
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Maslov's theory gBJ-gS]. A caustic arises due to singular- 
ities in the projection of some higher dimensional space 
to a lower one. In our case, the higher dimensional space 
is the phase space of null rays (with coordinates (x, p)) 
wherein the null rays follow trajectories with coordinates 
(x(<^), p(ip)) that do not intersect because their evolution 
is Hamiltonian [see (B7|]. 



Instead of forming the asymptotic series in l/u> for 
G(w;x, x'), Maslov's idea is to form the series for its 
Fourier transform with respect to at least one of its spa- 
tial coordinates. Then the Fourier transform operates 



as a canonical transformation in phase space. When 
the new momentum variables are projected out, caustics 
form at different locations than previously. This change 
of caustic location allows us, in principle, to construct a 
well-defined asymptotic series for the geometrical optics 
approximation of the Green function by using Maslov's 
theory whenever a caustic is about to be encountered 
(Fig. 23 ) . In this way, one may be able to derive di- 



rectly the phase shift through a caustic in (CI) without 
recourse to a matching argument. The application of 
Maslov's theory is outside the scope of this paper. 
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